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SUMMARY 

Supercavitating  bodies  can  achieve  very  high  speeds  under  water  by  virtue  of  reduced  drag:  with 
proper  design,  a cavitation  bubble  is  generated  at  the  nose  and  skin  friction  drag  is  drastically  reduced. 
Depending  on  the  type  of  supercavitating  vehicle  under  consideration,  the  overall  drag  coefficient  can 
be  an  order  of  magnitude  less  than  that  of  a fully  wetted  vehicle.  Slender-body  theory  and  boundary 
element  methods  are  two  modern  computational  methods  applied  to  the  design  of  supercavitating 
vehicles.  These  course  notes  present  recent  advances  in  the  theory  behind  these  two  computational 
approaches,  as  well  as  results  and  application  of  the  methods  to  the  simulation  and  control  of 
supercavitating  vehicles. 


1 INTRODUCTION 

Naval  hydrodynamics  and  the  marine  engineering  sciences  have  been  characterized  by  centuries  of 
evolutionary  innovation  punctuated  by  several  revolutionary  performance  improvements.  The 
relatively  continuous  improvement  in  early  paddle  wheel  technology,  for  example,  became 
obsolete  very  quickly  upon  introduction  of  the  screw  propeller  in  the  late  nineteenth  century.  The 
development  of  undersea  vehicles  has  seen  several  similar  paradigm  shifts,  such  as  the  introduction 
of  piston  engines  appropriately  packaged  to  fit  within  a typical  vehicle  envelope,  and  the 
application  of  acoustical  homing  techniques.  However,  until  recently,  there  has  been  little  interest 
in  very  high  speeds  for  undersea  applications. 

There  are  two  primary  impediments  to  increasing  underwater  vehicle  speeds:  current 
propulsion  capability  and  the  relatively  high  drag  of  current  underwater  vehicle  configurations. 
For  most  of  the  last  century  in  this  country,  relatively  little  research  has  been  devoted  to 
significantly  opening  the  performance  envelope  in  these  areas.  Over  the  last  decade,  however, 
high-speed  capability  has  received  increased  attention,  and  dramatic  advancements  have  been 
made.  These  achievements  have  been  supplemented  with  information  concerning  international 
development  (LEGI,  2000). 

The  history  of  hydrodynamics  research  displays  an  emphasis  on  eliminating  cavitation,  chiefly 
because  of  the  erosion,  vibration,  and  acoustical  signatures  that  often  accompany  the  effect.  The 
drag-reducing  benefits  of  cavitation,  however,  were  noted  during  the  first  half  of  the  last  century, 
and  have  received  significant  attention  over  the  last  decade. 

Supercavitation  is  a hydrodynamic  process  by  which  a submerged  body  is  almost  entirely 
enveloped  in  a layer  of  gas.  Because  the  density  and  viscosity  of  the  gas  is  dramatically  lower  than 
that  of  seawater,  skin  friction  drag  can  be  reduced  dramatically.  If  the  body  is  shaped  properly,  the 
attendant  pressure  drag  can  be  maintained  at  a very  low  value,  so  that  the  overall  body  drag  is  also 
reduced  significantly.  The  process  of  designing  a supercavitating  body  for  minimum  pressure  drag, 
while  addressing  issues  of  control  and  maneuvering,  is  greatly  aided  by  the  use  of  modem 
computational  methods. 

This  document  presents  details  of  two  computational  methods  that  have  been  applied 
extensively  during  the  past  decade  to  the  analysis  and  design  of  supercavitating  bodies:  boundary 
element  methods  and  slender-  body  theory.  Section  2 lays  the  groundwork  with  a discussion  of 
several  issues  pertinent  to  the  numerical  modeling,  including  a motivation  for  the  application  of 
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potential  flow  techniques  to  prediction  of  supercavitating  flows.  Section  3 presents  results  from 
several  applications  of  boundary  element  methods,  including  steady  axisymmetric  flows,  unsteady 
flows,  and  cavitating  fins.  Section  4 contains  a discussion  of  slender-body  theory  and  some 
pertinent  results.  In  section  5,  we  present  the  results  of  numerical  flight  simulation  for 
supercavitating  vehicles  that  incorporates  some  of  the  computational  methods  discussed  in  previous 
sections  in  determining  the  forces  experienced  by  the  vehicle  in  maneuvers. 

These  lecture  notes  were  derived  from  several  previous  publications.  The  nomenclature 
schemes  of  the  original  publications  have  been  retained,  but  each  sub-section  is  reasonably  self- 
consistent. 

2 PRELIMINARY  REMARKS 

2.1  Vehicle  Control 

A properly  shaped  supercavitating  body  will  have  a very  small  pressure  drag,  so  that  the  overall 
body  drag  is  also  reduced  dramatically  over  that  of  a fully  wetted  body.  However,  because  the 
center  of  pressure  is  typically  located  well  forward  with  respect  to  the  center  of  gravity,  control  and 
maneuvering  present  special  challenges.  Also,  whereas  a fully  wetted  vehicle  develops  substantial 
lift  in  a turn  due  to  vortex  shedding  off  the  hull,  a supercavitating  vehicle  does  not  develop 
significant  lift  over  the  gas-enveloped  surfaces.  This  requires  a different  approach  to  effecting 
hydrodynamic  control,  and  increases  the  relative  advantage  of  a banked  maneuvering  strategy. 

A supercavity  can  be  maintained  in  one  of  two  ways:  (1)  by  achieving  such  a high  speed  that 
the  water  vaporizes  near  the  nose  of  the  body,  generating  a cavity  that  grows  to  exceed  the  length 
of  the  body;  or,  (2)  by  supplying  gas  to  the  cavity  at  nearly  ambient  pressure  at  more  moderate  (but 
still  very  high)  speeds.  The  first  technique  is  known  as  vaporous  cavitation;  the  second  is  termed 
ventilation,  or  artificial  cavitation.  Schematic  views  of  two  types  of  notional  supercavitating 
bodies  are  shown  in  Figure  1 . Note  that  each  concept  involves  a cavitator  (in  this  case  a disk)  with 
a salient  edge  that  ensures  clean  cavity  formation  near  the  nose  of  the  body.  The  relatively  small 
diameter  of  the  cavitator  with  respect  to  that  of  the  vehicle  is  also  important:  it  is  this  feature  that 
allows  the  pressure  drag  to  be  maintained  at  a manageable  level.  An  important  step  in  the  design  of 
a supercavitating  body  is  selection  of  a cavitator  that  is  appropriately  sized  for  the  vehicle  and 
speed  of  interest. 

Control  of  supercavitating  vehicles  presents  special  challenges  not  normally  associated  with 
fully  wetted  vehicle  dynamics.  These  arise  from  the  absence  of  certain  physical  effects  such  as  lift 
on  the  body  and  from  the  presence  of  other  effects  such  as  the  nonlinear  interaction  of  the  control 
surfaces  and  the  body  with  the  cavity  wall. 

In  contrast  with  the  controlled  vehicle  depicted  in  Figure  lb,  a projectile  of  the  type  shown  in 
Figure  la  need  not  have  fins,  since  the  global  stability  of  the  trajectory  is  maintained  via  occasional 
tail-slap  contact  of  the  afterbody  with  the  cavity  boundary  (or,  at  extremely  high  speeds,  by  the 
forces  due  to  the  relative  motion  between  the  vapor  and  the  body).  Since  the  vehicle  shown  in 
Figure  lb  would  incorporate  some  sort  of  guidance  system,  fins  have  been  indicated  that  would 
provide  hydrodynamic  control  in  concert  with  actuation  of  the  cavitator.  The  blast  tube  (that  is,  the 
small-diameter  section  of  the  vehicle  aft  of  the  fins)  shown  in  Figure  lb  is  expected  to  represent  a 
more  optimal  shape  than  a cylinder  of  constant  radius  for  certain  classes  of  supercavitating 
vehicles:  if  range  is  critical,  then  drag  must  be  minimized,  which  (in  turn)  entails  minimization  of 
the  cavity  diameter.  The  step  at  the  after  end  of  the  mid-body  allows  for  this. 

2.2  Drag  Components  for  Submerged  Bodies 

For  simplicity,  consider  an  undersea  vehicle  with  constant  mass.  This  excludes  rocket  propulsion, 
but  is  illustrative  of  the  physics  of  supercavitation  drag  reduction.  The  following  analysis  of  the 
drag  force  is  equally  applicable  to  rocket-propelled  vehicles;  however,  the  basic  equation  of  vehicle 
motion  (equation  2.1)  must  be  modified  to  account  for  the  ejection  of  mass  (Greenwood,  1965). 
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Newton’s  second  law  of  motion  applied  to  undersea  vehicles  can  be  written  (Greenwood, 
1965) as 

T — D — — (mU) 

dt  (2.1) 

where  m is  the  vehicle  mass  plus  the  hydrodynamic  added  mass  (Newman,  1980)  and  U is  the 
velocity  of  the  vehicle  in  straight  and  level  flight  along  a trajectory  parallel  to  the  thrust  vector. 
Under  steady  conditions,  the  thrust  provided  by  the  propulsion  system,  T,  balances  the  vehicle 
drag,  D.  For  this  simple  case,  the  dynamics  and  hydrodynamics  of  the  vehicle  are  usually  analyzed 
in  a body-fixed  frame  of  reference,  as  shown  in  Figure  2. 

The  drag  represents  an  integration  over  the  vehicle  surface  of  all  components  of  stress  acting 
opposite  to  the  direction  of  travel  (Newman,  1980): 

D = O (pn  + n x + n x ) dS 

J7  V x y V z zx  ) 


where  p is  the  pressure,  xyx  and  x/x  are  the  components  of  the  shear  stress  in  the  x-direction,  and  n = 
(«x,  My,  mz)  is  the  unit  vector  normal  to  the  surface  of  the  vehicle. 

It  is  often  convenient  to  separate  the  total  drag  force  into  components  representing  the  pressure 
and  skin-friction  contributions: 

D=DP+DF  where  DP  is  the  pressure  drag  and  DF  is  the  skin-friction  drag. 

It  should  be  noted  that,  for  vehicles  traveling  close  to  a free  surface,  generation  of  gravity  and 
capillary  waves  is  associated  with  a modification  of  the  pressure  and  velocity  fields  on  the  vehicle 
hull,  causing  a drag  increase.  Also,  at  speeds  greater  than  the  lower  critical  Mach  number,  there  is 
a drag  increase  associated  with  the  shock  wave  system  (Ashley  and  Landahl,  1965).  In  most  cases 
of  interest,  the  primary  drag  component  affected  is  the  pressure  drag.  For  purposes  of  analysis, 
each  contribution  is  often  treated  as  a separate  component,  although  some  coupling  of  the  terms 
occurs.  Free-surface  drag  will  not  be  considered  in  the  following  discussion.  A compressible  flow 
formulation  using  slender-body  theory  is  presented  in  a later  section,  but  discussion  of  the 
associated  effects  on  drag  is  relegated  to  the  originating  publication,  Varghese,  et  al,  (1997). 

Denoting  the  fluid  density  as  p and  the  free  stream  speed  as  Ux,  the  drag  is  typically 
normalized  on  the  product  of  the  dynamic  pressure  of  the  free  stream,  /2  p ut> , and  a characteristic 
area,  A.  For  axisymmetric  undersea  vehicles,  A is  usually  chosen  as  the  maximum  sectional  area. 
Thus, 

C = ° 

° y2PUlA.  (2.3) 


2.3  Governing  Equations 

Neglecting  the  effect  of  gravity  (assuming  that  free-surface  and  stratification  effects  are 
unimportant  for  the  case  under  consideration),  the  equations  of  motion  of  a fluid  in  a volume  V in 
the  incompressible  limit  describe  conservation  of  mass: 

v • U = 0 (2.4) 


and  conservation  of  momentum: 


dt 


+ (U-V) 


U = -^-  + W2U, 


(2.5) 


where  V is  the  fluid  kinematic  viscosity.  The  momentum  conservation  equations  are  referred  to  as 
the  Navier-Stokes  equations. 
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The  relative  importance  of  the  terms  on  the  right-hand  side  of  equation  2.6  can  be  gleaned  by 
making  the  entire  system  dimensionless,  as  follows: 


~A  U ~AUj 
U= — t=-ss- 
U L 


„ A 

V = LV 


where  L is  some  length  characterizing  the  flow.  For  axisymmetric  undersea  vehicles,  this  length  is 
usually  chosen  as  the  maximum  diameter.  Making  these  substitutions  into  equations  2.5  and  2.6, 
rearranging,  and  dropping  the  tildes  gives 


V • V = 0 


(2.6) 


(i/.V) 


U=--VCv+  — V2U, 
2 p R 


where 


Cr. 


A P~  Poc 

~y2pul 


A 

— pressure  coefficient 


and 


(2.7) 


(2.8) 


J?  - . Reynolds  number, 

v 


(2.9) 


Prandtl  noted  that,  as  the  Reynolds  number  increases,  viscous  effects  become  increasingly 
confined  to  a thin  boundary  layer  near  the  body  surface  (Schlichting,  1979).  This  fact,  in 
conjunction  with  certain  mathematically  provable  laws  of  fluid  motion  (Yih,  1979),  allows  use  of 
potential  flow  techniques  in  predicting  the  motion  of  high-Reynolds-number  flows  past  bodies 
(Newman,  1980).  Under  these  conditions,  the  velocity  field  can  be  determined  based  solely  on 
kinematic  considerations  and  knowledge  of  conditions  at  the  boundaries  of  the  fluid.  Once  the 
irrotational  velocity  field  is  known,  an  auxiliary  equation  for  pressure  is  given  by  a simple  form  of 
Bernoulli’s  equation  (Aris,  1962): 


P 


(2.10) 


The  far  field  pressure,  p^,  is  simply  the  hydrostatic  pressure  at  the  operating  depth  of  the  vehicle. 

When  the  pressure  of  a liquid  is  reduced  below  some  threshold  level  at  constant  temperature  it 
vaporizes.  This  change  of  phase  is  known  as  cavitation  (Newman,  1980).  The  threshold  value  for 
a pure  liquid  undergoing  quasi-static  pressure  reduction  is  known  as  the  vapor  pressure,  pv. 
Cavitation  inception  is  often  enhanced  in  the  field  due  to  various  perturbations  of  ideal  conditions, 
including  the  amount  of  gases  absorbed  in  the  liquid,  contamination  by  other  liquids  or  solids,  or 
the  presence  of  bubbles.  The  vapor  pressure  provides  an  excellent  measure  of  the  cavitation 
inception  pressure  under  most  conditions  of  interest,  however. 

As  fluid  flows  past  a fully  wetted  body,  the  pressure  increases  to  stagnation  pressure  at  the 
forward  stagnation  point,  decreases  below  ambient  pressure,  then  increases  toward  stagnation 
pressure  over  the  afterbody.  The  actual  pressure  recovery  depends  on  the  development  of  the 
boundary  layer  and  the  extent  of  separation. 

If  the  pressure  falls  below  the  fluid  vapor  pressure  at  any  point,  a vapor  cavity  will  form.  The 
density  of  the  vapor  within  the  cavity  is  orders  of  magnitude  less  than  that  of  the  ambient  liquid. 
The  dimensionless  quantity  characterizing  the  tendency  of  a given  flow  to  cavitate  is  known  as  the 
cavitation  number: 

Poo-  Pc 

y2pui 


(2.11) 
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where  pc  is  the  cavity  pressure.  For  non-ventilated  flows,  the  cavity  pressure  is  approximately 
equal  to  the  vapor  pressure  of  the  ambient  liquid  under  most  conditions  of  current  interest. 

As  the  velocity  increases,  the  cavitation  number  and  the  minimum  pressure  on  the  body  surface 
decrease,  and  the  fluid  tends  to  cavitate.  Below  a certain  cavitation  number,  an  experimentally 
observed  fact  (one  that  may  be  explained  by  considering  the  various  terms  in  the  momentum 
equation,  equation  2.5)  is  that  the  velocity  within  the  vapor  cavity  is  very  small.  Consequently,  to 
lowest  order,  the  pressure  within  the  vapor  cavity  is  constant,  with  a value  nearly  equal  to  the  vapor 
pressure  of  the  fluid.  At  high  Reynolds  number,  the  shear  layer  forming  the  cavity  boundary  is 
thin.  In  that  case,  potential  flow  techniques  may  be  used  to  predict  such  flows  (Newman,  1980).  It 
should  be  noted  that  recent  results  of  high-Mach-number  research  indicate  that  the  assumption  of 
negligible  velocities  within  the  cavity  is  not  valid  at  very  high  speeds  (say,  at  Mach  numbers 
approaching  unity).  However,  for  cavity-riding  vehicles,  the  pressure  can  be  taken  as  constant 
within  the  vapor  cavity. 

In  terms  of  velocity,  the  boundary-value  problem  describing  the  ideal  superca  vita  ting  flow 
shown  in  Figure  3 is  given  by  the  field  equation  describing  conservation  of  mass  (equation  2.4) 
along  with  boundary  conditions  on  velocity: 

U-n  = 0 on  Sw  u Sc  (2  i2i 


U-s  = U^-yl  1 + o on  ,SC 


(2.13) 


where  Sw  is  the  wetted  surface  of  the  body  and  Sc  is  the  cavity  boundary. 

The  dynamic  condition  results  from  applying  Bernoulli’s  equation  (equation  2.10)  along  the 
cavity  streamline,  under  the  constraint  that  the  pressure  within  the  cavity  is  constant  with  a value 
equal  to  the  cavity  pressure,  pc.  Various  techniques  exist  for  solving  such  a problem,  including 
panel  methods  and  slender-body  theory  (Kuethe  and  Chow,  1976).  These  two  approaches  will  be 
discussed  in  the  following  sections. 


3 APPLICATION  OF  THE  BOUNDARY-ELEMENT  METHOD 

This  section  describes  application  of  the  boundary  element  method  to  prediction  of  supercavitating 
flows.  Three  sample  problems  are  discussed:  (1)  steady,  subsonic  flow  around  an  axisymmetric 
disk;  (2)  unsteady,  subsonic  flow  around  general  cavitator  shapes;  and  (3)  steady  flow  around 
cavitating  control  effector  fins  for  application  to  supercavitating  vehicle  flight  control. 


3.1  Steady  Axisymmetric  Flows 

This  sub-section  is  derived  from  Kirschner,  et  al,  (1995). 

The  physical  problem  is  shown  in  the  schematic  drawing  of  Figure  4.  The  body  is 
axisymmetric  and  oriented  parallel  to  the  flow.  The  geometries  under  consideration  in  these  lecture 
notes  include  only  those  for  which  a salient  cavity  detachment  point  can  be  identified,  although  the 
solution  techniques  described  could  be  extended  to  more  general  shapes. 

The  position  vector  in  a meridonal  plane  is  denoted  as  x = xi  +rj . The  inner  boundary  of  the 
flow  problem  is  the  union  of  the  wetted  portion  of  the  body  and  the  boundary  of  the  cavity: 

S = SBuSc.  (3J) 

The  lengths  of  the  body  and  cavity  are  denoted  as  fB  and  /c , respectively.  The  total  length  of  the 
combined  body-cavity  system  is  denoted  as  A 

Only  steady  flows  are  considered  in  this  sub-section.  The  difficulties  of  representing  cavity 
closure  under  this  assumption  are  discussed  below. 

Unit  vectors  normal  and  tangent  to  the  body-cavity  surface  in  a meridonal  plane  are  denoted  as 
n = nxi  + nrj  and  s = sj  + sT  j , respectively.  The  unit  normal  vector  is  directed  into  the  liquid;  the 
unit  tangent  vector  is  directed  positive  aft.  Um  = Umi  is  the  free  stream  velocity.  The  total  velocity 
at  any  point  outside  of  the  body-cavity  surface  is  decomposed  as 
U = UeJ  + u 


(3.2) 
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where  u(x)=ux  i+urj  is  the  disturbance  velocity.  As  discussed  above,  to  lowest  order  the  cavity 
pressure,  pc 

Potential  Flow  Formulation 

The  fluid  is  assumed  to  be  incompressible  and  the  flow  is  taken  to  be  irrotational.  The  last 
assumption  guarantees  the  existence  of  a velocity  potential.  Under  this  condition,  the  flow  field  is 
governed  by  Laplace's  equation, 

V20  = </>„  + 0rr  +^~  = 0 

r , (3.3) 

where  d>  is  the  total  potential. 

A disturbance  potential,  <p , can  be  defined  from  the  total  potential  by 
0=Uao-X+<!>  (3  A\ 


so  that  the  disturbance  velocity,  u,  is  given  by  the  gradient  of  the  disturbance  potential.  The 
disturbance  potential  also  obeys  the  Laplace  equation.  This  boundary  value  problem  is  shown  in 
the  schematic  drawing  of  Figure  5. 


Kinematic  Condition 

The  requirement  that  the  flow  be  tangential  to  the  body  and  cavity  surfaces,  equation  2.12, 
combined  with  the  definition  of  the  potential  in  equation  3.4,  gives  the  following  kinematic 
condition  on  the  body-cavity  surface,  S : 


d(p 

dn 


—n  U^onS 


(3.5) 


Dynamic  Condition  On  The  Cavity 

Bernoulli's  equation  applied  along  the  cavity  surface  can  be  expressed  as  follows: 
Cp  = -a  on  Sq 


The  magnitude  of  the  total  velocity  on  the  cavity  surface  is  given  by  equation  2.13.  From  this,  the 
following  expression  for  the  disturbance  potential  may  be  applied  on  the  cavity  surface: 

= LL Vl+  c - on  Sc 

ds  . (3.7) 


Cavity  Termination 

For  a given  wetted  body  shape  and  cavity  length,  Laplace’s  equation  for  the  disturbance  potential 
can  be  solved  subject  to  the  kinematic  condition,  equation  3.5,  and  the  dynamic  condition, 
equation  3.7.  A cavity  termination  model  is  also  required. 

In  this  work,  a modified  Riabouchinsky  cavity  termination  model  is  used.  The  cavity  is  closed 
with  a Riabouchinsky  wall  or  endplate  extending  from  the  axis  of  symmetry  to  the  point  at  which 
the  cavity  streamline  becomes  perpendicular  to  the  axis.  Along  the  endplate,  denoted  as  Sw , the 
kinematic  condition  is  satisfied.  In  the  current  case,  since  the  flow  is  axisymmetric,  the  endplate  is 
always  perpendicular  to  the  axis  of  symmetry. 

It  should  be  noted  that  this  termination  model  has  been  chosen  as  a compromise  between 
simplicity  and  physical  validity.  The  reentrant  jet  model  is  a closer  representation  of  the  physical 
flow,  but  is  more  difficult  to  implement  numerically.  The  reentrant  jet  model  is  discussed  in 
further  detail  below. 

Solution  Procedure  using  Boundary  Elements 

The  disturbance  velocity  potential  in  equation  3.4  satisfies  Green's  third  identity,  a Fredholm 
integral  equation  of  the  second  kind,  along  the  body,  cavity,  and  endplate  (Lamb,  1945).  Thus,  for 


any  field  point  x on  the  body-cavity  surface,  the  disturbance  velocity  potential  can  be  computed 
from 
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2tt0(x)  = -O  4>(x')-^G(x;x')—^-4>(x')G(x;x') 
JJ  an  an 


dS(x') 


(3.8) 


where  the  integration  is  performed  over  source  points,  x\  along  the  entire  body,  cavity,  and 
endplate  surface.  Green's  function  for  this  problem,  GCxjx'h  is 


G(x;x') 


(3.9) 


Along  the  body  and  the  endplate  surfaces,  the  source  strength  distribution,  d(j>/dn,  is  known 
and  given  by  equation  3.5.  Along  the  cavity  surface,  an  expression  for  the  dipole  strength 
distribution,  0 , can  be  found  by  integrating  equation  3.7: 

0(i’)  = 0(i’o)  + f/„-\/l  + (T  (s-Sq)-U^  (jC-JCq)  on  Sc- 

(3.10) 


where  s0  is  the  arc  length  coordinate  at  an  arbitrary  fixed  point  on  the  body-cavity  surface,  S . 

To  solve  for  the  dipole  strength  distribution,  0 , along  the  body  and  endplate  surfaces,  and  the 
source  strength  distribution,  d$/dn,  along  the  cavity  boundary,  the  computational  domain  is 
discretized  into  MB  panels  along  the  body,  Mc  panels  along  the  cavity,  and  Afw  panels  along  the 
endplate.  Unequal  spacing  of  the  panels  is  utilized  to  ensure  that  more  panels  are  concentrated 
where  large  velocity  gradients  are  expected. 

The  source  and  dipole  strengths  are  taken  to  be  constant  on  each  panel.  Integral  equation  3.8  is 
satisfied  at  control  points  on  every  panel.  These  are  chosen  such  that  each  trapezoidal  panel  is 
divided  into  equal  areas  upstream  and  downstream  of  its  control  point. 

The  unknowns  to  be  determined  are  therefore 


• Mb  dipole  strengths  on  the  body  panels; 

• Mq  source  strengths  on  the  cavity  panels; 


• Mw  dipole  strengths  on  the  endplate  panels;  and, 

• the  dimensionless  cavity  surface  velocity,  Vl  + 


There  are  MB+MC  + Mw  + 1 quantities  to  be  determined, 
equations.  The  zero  net  flux  condition. 


dS  = 0 


Equation  3.8  gives  MB+MC+MW 


(3.11) 


is  used  as  another  condition  to  form  a determinate  system. 

Cavity  Shape  For  Fixed  Cavity  Length 

For  the  body  shapes  considered  in  these  lecture  notes,  the  cavity  detachment  point  is  known.  For  a 
specified  cavity  length,  a cavity  profile  is  assumed,  and  the  potential-based  boundary  element 
method  used  to  solve  for  the  dimensionless  cavity  surface  velocity,  Vl  +tr , and  the  source  and 
dipole  strengths.  A new  estimate  of  the  cavity  profile  is  then  determined  by  applying  the  kinematic 
condition  on  the  cavity  surface,  5c  . The  boundary  value  problem  is  then  re-solved,  and  the  process 
repeated  until  the  cavity  profile  converges. 

The  kinematic  condition,  equation  3.5,  can  be  written  as 
nx  ( mx  +{/«,)=  — nTuT . 


(3.12) 
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For  example,  if  the  body  cavity  surface  profile  is  defined  by  y = b(x),  then 

db  «v  iir 

dx  nr  ux  + tfy 


(3.13) 


along  the  cavity  surface. 

Prior  to  convergence  of  the  cavity  profile,  the  kinematic  condition  on  the  cavity  boundary  has  not 
been  imposed.  Thus,  the  object  is  to  change  the  cavity  shape  such  that  the  flux  across  the  boundary 
is  minimized.  Two  methods  have  been  tested  for  modifying  the  cavity  shape  including  (1) 
modification  of  the  profile  offset;  and,  (2)  sequential  modification  of  each  panel  angle  followed  by 
re-scaling  of  the  cavity  to  its  original  length.  In  practice,  the  latter  approach  has  proven  to  be  more 
general  and  robust  than  the  former. 

Pressure  and  Drag 

Once  the  potential  flow  problem  outlined  in  the  previous  section  has  been  solved,  the  pressure 
coefficient  on  the  body  surface,  SB,  can  be  computed  as 


Cp=l- 


'U^ 

U 


= 1- 


Sx  + " 


U 


d(j) 

ds 


(3.14) 


The  body  drag  is  the  integral  of  the  x -component  of  the  differential  pressure  force  over  the 
wetted  portion  of  the  body  and  the  base: 


D= 


II 


Pr>x 


dS  + kB 2 


Pc 


(3.15) 


Where  B is  the  radius  of  the  cavitator  at  the  cavity  detachment  point.  The  drag  coefficient  is 
computed  as: 


CD  = cr  - 


Sb 


«VC„ 


dS 


(3.16) 


Selected  Results  - Steady  Flows 

Sample  results  of  the  fully  axisymmetric  boundary  element  method  for  incompressible  flow  are 
presented  in  Figure  6 through  Figure  8.  For  additional  results  and  discussion,  refer  to  Kirschner,  et 
al,  (1995). 

The  body-cavity  profiles  predicted  for  a 15°  cone  at  several  cavitation  numbers  are  shown  in 
Figure  6.  The  predicted  drag  coefficients  at  each  cavitation  number  are  plotted  in  Figure  7 and 
compared  with  the  empirical  curve  fit  of  May  (1975): 

CD  = 0.151  +0.587(7  (3  17) 

It  can  be  seen  that  the  trend  is  very  accurate.  Since  the  publication  of  Kirschner,  et  al  (1995),  a 
coding  error  was  discovered  which,  when  corrected,  virtually  eliminated  the  offset  seen  in  Figure  7. 

The  predicted  maximum  cavity  diameter  is  compared  with  the  empirical  curve  fit  of  May 
(1975)  in  Figure  8.  The  following  equation  represents  the  curve  fit: 

— = 1. 28cF0'91 

CdI  b ) (3.18) 

The  left-hand  side  of  equation  3.18  is  the  function  plotted  on  the  ordinate  ofFigure8.  The 
prediction  appears  to  be  excellent.  Note  from  the  figures  that,  if  the  predicted  drag  coefficient  used 
in  equation  3.18  is  replaced  with  that  given  by  equation  3.17,  the  prediction  of  values 
characterizing  the  cavity  geometry  is  improved.  It  should  be  noted  that  May’s  equation  over- 
predicts experimental  cavity  lengths  at  higher  cavitation  numbers. 
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Reentrant  Jet  Termination  Model 

A termination  model  that  better  represents  the  physics  for  typical  supercavitating  flows  of  interest 
is  the  reentrant  jet  model,  a description  of  which  may  be  found  in  Tulin  (1964).  The  current 
formulation  is  based  on  Uhlman,  et  al,  (1998). 

For  axisymmetric  flows,  the  presence  of  the  reentrant  jet  may  be  handled  quite  simply  by 
truncating  an  axisymmetric  reentrant  surface  at  a jet  face,  and  assuming  that,  at  this  axial  location, 
this  surface  is  asymptotic  to  a cylinder  of  constant  diameter  and  the  cavity  surface  velocity  has 
approached  a constant  limit.  In  the  absence  of  gravity,  and  assuming  that  the  truncated  jet  face  is 
perpendicular  to  the  velocity,  the  conditions  at  the  jet  face  are 

30  TT  r 

— = f700Vl  + cx 

3«  (3.20) 

and 

3<t> 

— = 0. 

ds  (3.21) 


With  the  disturbance  potential  defined  by  equation  3.4,  the  conditions  to  be  applied  at  the  jet  face 
become 


— = uJ>fi+o+i) 

dn  v 1 


and 


ds 


(3.22) 


(3.23) 


Note  that  for  two-dimensional  reentrant  jets  with  lift  the  conditions  are  more  complicated.  Figure  9 
shows  a sample  computed  cavity  shape  for  a two-dimensional  flat  plate  cavitator  at  an  angle  of 
attack. 


3.2  General  Incompressible  Flows 

This  section,  which  is  based  on  Kring  et  al  (2000),  presents  some  initial  results  for  a fully  three- 
dimensional  time-domain  simulation  of  the  flow  about  general  supercavitating  bodies.  The  intent 
of  the  effort  is  to  study  various  ventilated  cavitators,  lifting  appendages,  and  bodies  that  may 
penetrate  the  cavity  wall.  The  time-domain  approach,  useful  for  the  formulation  of  the  nonlinear 
boundary  conditions  on  the  cavity  free  surface,  will  also  allow  the  application  of  this  method  to 
flight  simulation  and  control  studies  (see  section  5). 

In  this  fully  three-dimensional  initial  value  problem  for  a moving,  cavitating  body,  several 
boundary  conditions  must  be  specified.  These  include  the  body  boundary  condition,  nonlinear 
cavity  conditions  on  the  moving  free  surface,  a reentrant  jet  closure,  and  a shed  wake  sheet.  At 
each  time  step,  these  conditions  are  integrated  both  in  time  and  in  space  and  combined  with  the 
governing  Laplace  equation  to  produce  a boundary  value  problem  that  is  formulated  as  a boundary 
integral  equation  of  the  mixed  Neumann-Dirichlet  type. 

The  numerical  scheme  in  this  paper  has  been  adapted  from  a boundary-element  method 
developed  for  ffee-surface  gravity  waves  interacting  with  submerged  and  floating  bodies  (Kring,  et 
al,  1999),  with  guidance  for  the  cavity  boundary  and  jet  closure  conditions  from  Uhlman  (1987  and 
1989)  and  Fine  and  Kinnas  (1993).  The  free  surface  scheme  was  designed  in  accordance  with  the 
numerical  stability  analysis  of  Vada  and  Nakos  (1993).  A key  result  of  that  analysis  was  the 
demonstration  that  both  high-order  elements  and  mixed  explicit-implicit  time  marching  are  needed 
for  numerical  stability.  The  boundary-element  method  employed  here  is  a high-order,  B-spline- 
based  approach  developed  originally  by  Maniar  (1995)  with  a so-called  fast  or  accelerated  solver 
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developed  by  Phillips  and  White  (1997)  and  applied  in  Kring,  et  al,  (1999b).  Whereas  direct 
solution  of  full  matrices  requires  computational  effort  proportional  to  the  cube  of  the  number  of 
unknowns,  N,  — that  is,  effort  of  order  N 3,  or  0(jV3)  - and  other  traditional  methods  may  show 
improved  performance  of  0(iV2),  the  accelerated  method  discussed  in  this  section  is  0(N). 

Figure  10  illustrates  the  potential  flow,  initial  boundary-value  problem  for  a cavitating  disk 
with  a trailing  cavity,  reentrant  jet  closure,  and  shed  wake.  The  body  may  have  any  shape  and  the 
method  may  be  extended  to  add  lifting  appendages  and  planing  bodies. 

In  order  to  provide  a computational  window  to  track  this  moving  body,  the  three-dimensional 
problem  is  defined  with  respect  to  a mean-body  frame  of  reference,  x . Points  in  the  inertial  frame, 
x, , are  related  to  the  mean-body  frame  by  the  relation: 


x,  =x  + jU(x,T)dT, 

k) 


(3.24) 


where  t is  time.  The  motion  of  the  mean-body  frame,  U(x,t) , is  a function  of  space  since  it  may 
be  rotating  with  respect  to  the  inertial  frame,  so  that 


U (x,  t)  = U (0,  t ) + Q(t  )x, 


(3.25) 


where  Q (t)  is  a rotational  transformation  matrix  and  U(0,t)  is  the  translation  at  the  origin.  The 
mean-body  frame  is  defined  about  the  assumed  mean  position  of  the  disk  cavitator,  so  that  the  disk 
tends  to  oscillate  around  x = 0 . 

Initial  Conditions 

The  problem  is  started  impulsively  from  an  assumed  initial  geometry  and  distribution  of  potential. 
Two  initial  conditions  have  been  considered.  In  the  first  case,  the  potential  on  the  cavity  is  set  to 
zero  on  an  assumed  surface.  Physically,  this  can  be  considered  as  an  artificial  pressure  and  cavity 
location  that  is  released  at  t = 0.  If  the  cavitator  disk  is  set  in  steady  motion,  the  transients  due  to 
this  artificial  start-up  should  eventually  decay.  The  second  case  involves  solution  of  a Neumann 
problem  in  which  the  cavity  is  considered  to  be  a fixed  wall.  The  potential  from  this  solution  is 
then  used  as  the  initial  condition  for  the  time-domain  simulation.  The  intent  is  to  reduce  the  size  of 
the  startup  transient. 

Governing  Equation 

The  potential  flow  exterior  to  the  surface  S (the  union  of  the  cavity  S > and  S ; the  body,  SB ; 
and  the  jet,  Sj ) is  governed  by  the  Laplace  equation  with  allowance  for  the  wake  cuts  collapsed 
onto  the  surface  Sw  . Hence, 


V20  = O 


(3.26) 


where  0 is  the  disturbance  potential. 

Green’s  third  identity,  equation  3.8,  is  applied  in  a manner  similar  to  that  described  for  steady 
axisymmetric  flows.  In  this  case,  however,  we  explicitly  include  an  integration  over  the  wake 
surface,  Sw,  as  follows: 


2n(j){x)  4-  jj  (0(x  ')G(x;x ')  -0„  (x  ')G„  (x;x  '))dx'  + Jj  A0(x  ')G„(x;  x j dx ' = 0 


(3.27) 


where  the  Green’s  function  is  defined  by  equation  3.9. 
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Cavity  Conditions 

The  cavity  conditions  on  £ + and  S are  posed  in  the  curvilinear  frame  of  reference  defined  by 

the  parametric  surface,  £,  (u,  v,t) . The  parametric  coordinates,  (u,v),  are  the  presumed 

streamwise  and  circumferential  directions,  respectively.  This  semi-Eulerian  formulation  is 
important  for  a numerically  consistent  representation  of  the  boundary  condition  on  the  moving 
surface.  Surface  potentials,  (p(u,v,t)  and  (p0(v,t) , are  defined  through  the  decomposition  of  the 
perturbation  potential: 

The  first  part  of  the  surface  potential,  (p , satisfies  the  dynamic  cavity  condition  implicitly  in  the 
numerical  scheme.  The  second  part,  (pa,  which  is  constant  along  the  presumed  streamline 

direction,  u , satisfies  continuity  of  the  disturbance  potential  at  the  cavity  detachment  locus  and 
appears  explicitly  in  the  dynamic  cavity  condition  on  the  leading  cavity  surface. 

Applying  the  chain  rule  to  the  decomposition  leads  to  relations  for  the  time-derivative  and 
gradient: 

d(j)  _d((p  + (pQ) 
dt  dt 


-v,  (<?+«,)■  v„f 


and 


V*  = Vk(?>  + *.)  + fpi 


(3.29) 


(3.30) 


where  h is  the  unit  vector  normal  to  the  surface  and  Vw  is  the  surface  gradient  operator  defined  by 


V,/ 


OT  r a / a a\  a~i  CJT  r~  a /a  a\  a-] 

— u-(u-v)v  +~  v-lu-v)u 
_ dl  l L J ov  L J 


MXV 


(3.31) 


Here  u and  v are  unit  vectors  that  are  non-orthogonal  in  general. 


Kinematic  Boundary  Condition 

The  total  time  derivative  of  the  function  that  defines  the  cavity  surface  must  vanish.  Applying  the 
Gallilean  transform  from  the  inertial  to  the  mean-body  frame  and  defining  the  normal  coordinate, 

n , and  the  surface  displacement  in  the  normal  direction,  h(u,v,t)  = %(u,v,t)-n(u,v,t ) , yields  the 
kinematic  condition: 


^._u-vuh=-vu((p+(p0yvuh-u-h+^- 

ot  an 


(3.32) 


Dynamic  Boundary  Condition 


The  dynamic  condition  specifies  a constant  pressure  on  the  cavity  surface  at  all  times.  The 
derivation  of  this  condition  begins  with  the  unsteady  Bernoulli  condition  in  the  inertial  frame 
(referencing  atmospheric  pressure  as  zero).  With  z = 0 representing  the  mean  ocean  surface: 


Pc=-P 


— + -V0-V0  + gz 
dt  2 Y V & 


J 


(3.33) 


where  pc  is  the  cavity  pressure. 
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Applying  the  Gallilean  transform  from  the  inertial  to  the  mean-body  frame  and  substituting 
relations  for  the  surface  potential  yields 


d(p 


U-Vu(p 


dt 


+ U-'Vu(p()  + 


(3.34) 


Cavitator  Wetted-Surfaee  Condition 

A Neumann  body-boundary  condition  similar  to  equation  3.5  is  applied  to  the  cavitator  wetted 
surface,  SB . 


Jet  Condition 

The  reentrant  jet  is  artificially  truncated  (cut  perpendicular  to  the  streamline)  and  assumed  to  have 
reached  a steady  diameter  with  no  cross-flow  components  and  a pressure  equal  to  the  cavity 
pressure  across  the  face  of  the  cut.  Thus,  it  has  an  imposed  flux  and  an  imposed  pressure. 
Although  this  is  intended  to  be  a time-accurate  formulation  for  the  cavity  problem,  a steady-state 
jet  condition  has  been  used  to  date. 

We  start  with  the  steady,  dynamic  boundary  condition:  equation  3.34  with  the  time  derivative 
set  to  zero.  If  there  is  no  cross-flow,  then  V (j)  —(j)nn  on  the  jet.  Therefore,  the  jet  condition 
becomes 


t = gz+£_ 

dn  2{dn  ) p , (3.35) 


where  Un  -U -n. 


Choosing  the  positive  root,  consistent  with  the  cavity  conditions,  the  jet  flux  then  becomes: 


P 


A 


(3.36) 


With  no  cross-flow,  the  potential  across  the  jet  face  must  be  constant  at  a given  depth  and  is  set 
equal  to  the  potential  on  the  adjacent  cavity. 


Cavity  Detachment  and  Kutta  Conditions 

The  additional  potential,  (p0 , on  the  leading  cavity  surface,  S + , that  detaches  from  the  cavitator 

body  boundary  is  specified  by  a detachment  condition  across  the  body-cavity  intersection,  x \BE  . 
With  the  further  constraint  that  the  first  surface  potential  vanishes  at  the  cavity  detachment  line, 
(p  (u  \BE  5 v,  t ) = 0 , the  additional  potential  is. 


% (v,t)  = $(x  | BE,t) 


(3.37) 


which  represents  the  cavity  detachment  condition.  A cavity  trailing  edge  is  defined  at  the  junction 
of  the  leading  and  lagging  cavity  surfaces  where  the  potential  jump  across  the  wake  is  specified  as 


A(f)  TE  (V,f)  {(p  \TE  (p  + o 


(3.38) 


which  comprises  the  Kutta  condition.  For  a time-dependent  wake,  a nonlinear  evolution  equation 
should  govern  the  propagation  of  the  shed  vorticity  downstream,  but,  to  date,  a simplifying 
assumption  has  been  used:  The  wake  position  is  fixed  and  a steady  state  asymptotic  distribution  is 
assumed  for  the  shed  vorticity. 
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Numerical  Implementation 

The  scheme  begins  with  initial  conditions  for  the  geometry  and  the  surface  potential,  (p , at  time 
t = 0.  The  detachment  conditions,  the  Kutta  conditions,  and  the  mixed  boundary-integral  equation 
are  then  solved  simultaneously  to  find  the  normal  velocity  on  the  cavity  surface  and  the  potential 
on  the  body.  This  is  referred  to  as  the  mixed  BVP,  which  will  be  solved  at  each  time  step  as  the 
surface  evolves.  The  cavity  conditions  and  the  mixed  BVP  must  be  integrated  in  time  and  space. 

After  the  impulsive  start,  each  step  in  the  temporal  integration  begins  by  first  satisfying  the 
kinematic  and  dynamic  boundary  conditions  to  update  the  geometry  and  surface  potential  on  the 
cavity.  The  mixed  BVP  is  then  solved  to  update  the  potential  on  the  body  and  the  normal  velocity 
on  the  cavity. 

To  summarize,  at  the  start-up  (7  = 0" ) : 

0.  Given  initial  conditions  <p°  ,<f 0 on  Sc  solve  mixed  BVP  for  0(1°  on  Sc  . 

At  each  subsequent  time-step: 

1 . Determine  body  motion:  impose  or  solve  equation  of  motion. 

2.  Integrate  the  kinematic  condition  on  the  cavity  to  find  "+1  on  Sc  . 

3.  Update  the  cavity  and  body  geometry. 

4.  Integrate  the  dynamic  condition  on  the  cavity  to  find  <p',+l  on  Sc . 

5.  Set  the  body  and  jet  conditions  for  (pj : on  SB  and  S7 . 

6.  Solve  mixed  BVP  to  find  (j)n"+' ,tpQ"  1 on  Sc  and  (pn+l  on  SB . 

The  continuous  problem  is  discretized  by  a set  of  high-order  B-spline  patches  that  represent  the 
geometry  of  the  body  and  cavity  and  another  set  that  represent  the  flow  quantities  of  potential,  (p , 

normal  velocity,  (j)n , and  cavity  displacement,  h . 

Third-  to  fifth-order  splines  and  Galerkin  schemes  are  used  to  spatially  integrate  both  the  cavity 
conditions  and  the  mixed  BVP.  The  use  of  a high-order  discretization  has  been  shown  to  be  stable 
for  the  free-surface  gravity  wave  problem  where  more  typical  flat  panel  methods  have  failed. 
Initial  work  with  low-order,  three-dimensional  panel  methods  indicated  that  this  would  hold  for 
cavity  free-surface  flows  also. 

The  kinematic  cavity  condition  is  integrated  in  time  using  a forward  (explicit)  Euler  method 
and  the  dynamic  cavity  condition  is  integrated  in  time  using  a backward  (implicit)  Euler  method. 
For  the  dynamic  condition,  only  the  linear  terms  for  the  surface  potential  are  treated  implicitly. 
This  mix  of  explicit  and  implicit  methods  was  selected  based  on  a stability  analysis  for  wave 
propagation  over  discretized  free  surfaces. 

The  geometry  and  surface  potential  are  represented  by  B-spline  basis  functions  in  space  that 
overlap  with  their  neighbors.  These  conditions  produce  a banded  system  of  linear  equations  that 
must  be  solved  at  each  time-step. 

The  unknown  normal  velocity  on  the  cavity  and  potential  on  the  body  are  obtained  from  the 
solution  of  the  boundary  integral  equation,  which  contains  the  detachment  and  Kutta  conditions. 
On  the  cavity,  the  disturbance  potential  at  any  instant  is  the  sum  of  the  known  surface  potential  and 
an  unknown  additional  potential. 

The  solution  of  the  integral  equation  at  each  time-step  represents  the  largest  computational 
burden  in  the  problem.  To  maintain  high  local  accuracy  and  satisfy  numerical  stability 
requirements,  a high-order,  B-spline,  Galerkin  boundary-element  method  is  used  for  this  integral 
equation  and  two  solvers  are  adopted.  The  first  is  a traditional  dense  matrix  solver  that  requires 
O (N2)  computational  memory  and  operations.  The  second  is  an  accelerated  solver,  using  a 
precorrected  fast  Fourier  transform  (FFT)  algorithm  that  is  similar  to  the  fast  multipole  method. 
This  involves  a sparse  matrix  requiring  only  0(A)  computational  memory  and  operations. 


9-14 


Selected  Results  - Genera!  Flows 

To  date,  the  high-order  boundary-element  method  has  been  adapted  for  the  mixed  BVP  with  the 
detachment  and  Kutta  conditions.  In  Figure  1 1 an  initial  solution  of  the  problem  at  the  startup  has 
been  computed.  In  this  figure,  the  high  negative  values  of  the  normal  velocity  at  the  trailing  end  of 
the  cavity  show  that  the  assumed  geometry  is  too  small  and  the  cavity  will  lengthen  along  the  side 
in  the  next  time  step.  Work  on  the  complete  simulation  is  in  progress. 

This  boundary-element  method  is  numerically  convergent,  has  been  compared  to  analytic 
solutions  for  simple  geometries,  and  has  been  validated  experimentally  in  its  original  ship 
seakeeping  manifestation,  but  it  is  too  early  to  judge  the  performance  of  the  cavity  time-marching 
scheme  and  the  wake  sheet. 

3.3  Cavitating  Fins 

In  this  section,  we  describe  a tool  used  to  predict  the  flow  around  cavitating  fins,  particularly  fins 
with  wedge-shaped  cross-sections.  The  tool  uses  a boundary  element  method  to  predict  six-degree- 
of-freedom  forces  for  steady  flight.  The  computed  forces  have  been  used  as  input  to  guidance  and 
control  simulations,  as  described  below.  The  basic  computational  approach  is  summarized  in  Fine 
and  Kinnas  (1993). 

A simple  wedge-shaped  fin  geometry  was  considered  for  this  investigation.  It  would  be  easily 
fabricated  and  appropriate  to  operation  in  the  supercavitating  regime,  but  is  probably  not  optimal, 
although  the  wedge  shape  provides  good  strength  characteristics.  It  is  also  convenient  for 
comparison  of  the  boundary-element  results  with  preliminary  predictions  based  on  two- 
dimensional  theory  such  as  that  of  Uhlman  (1987,  1989)  or  Wu  (1956),  or  the  various  experimental 
data  presented  in  May  (1975).  Also,  for  small  wedge  angles,  the  unsteady  effects  associated  with 
the  transition  between  partial  cavitation  and  supercavitation  are  confined  to  a very  limited  range  of 
deflection  angles. 

In  computing  the  forces,  a high  level  of  detail  was  required  over  even  the  limited  operational 
range  considered,  because  several  cavity  detachment  modes  must  be  taken  into  account.  The  forces 
acting  on  cavitating  fins  are  complicated  by  these  different  flow  regimes.  For  wedge-shaped  fins, 
four  regimes  must  be  considered.  If  the  flow  is  symmetric  about  the  symmetric  fin,  and  the 
cavitation  number  is  low  enough,  a cavity  develops  at  the  base  of  the  wedge.  Of  course  the  fin 
does  not  develop  lift  in  this  case.  At  small  angles  of  attack  and  moderately  low  cavitation  numbers, 
a partial  cavity  develops  at  the  leading  edge  separately  from  the  base  cavity.  At  larger  angles  of 
attack  or  lower  cavitation  numbers,  the  cavity  originating  at  the  leading  edge  becomes  large  enough 
to  merge  with  the  base  cavity,  forming  a supercavity  that  envelopes  all  the  surfaces  except  the 
pressure  face.  The  lift  curve  slope  in  the  supercavitating  regime  is  less  than  in  the  partially 
cavitating  regime.  The  final  regime  (which  has  not  been  modeled  for  the  current  effort)  involves 
very  high  cavitation  numbers.  In  that  case,  the  cavities  are  eliminated  and  replaced  with  separation 
regions.  The  three  regimes  are  shown  in  Figure  12. 

Brief  Description  of  LScav 

ETC’s  Lifting  Surface  Cavitation  code  (LScav)  computes  the  first  iteration  of  the  fully  nonlinear 
cavity  solution  for  a general  three-dimensional  wing  attached  to  a non-lifting  “center-body”  at 
arbitrary  cavitation  number.  The  flow  is  assumed  to  be  steady,  inviscid,  incompressible,  and 
irrotational  (outside  of  a thin  vortex  sheet  trailing  behind  the  lifting  surface).  The  solution  is  found 
by  applying  a disturbance-potential-based  panel  method  that  involves  a numerical  solution  of 
Green’s  third  identity.  The  method  is  described  in  the  context  of  cavitating  propellers  in  Fine 
(1992). 

The  wing  and  center-body  are  assumed  to  be  symmetric  about  a vertical  plane  intersecting  the 
body  centerline.  The  cavity  solution  represents  the  first  iteration  of  a fully  nonlinear  solution.  In 
order  to  understand  the  difference  between  the  first  iteration  solution  and  the  fully  nonlinear 
solution,  consider  a two-dimensional  hydrofoil  and  assume  that  the  cavity  length  is  known.  Since 
the  cavity  length  is  a single-valued  function  of  the  cavitation  number  for  all  stable  two-dimensional 
cavities,  the  cavitation  number  may  not  be  specified  in  this  canonical  problem  and,  in  fact,  must  be 
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considered  an  unknown  to  be  determined  as  part  of  the  solution.  To  solve  the  nonlinear  problem, 
one  must  first  guess  the  location  of  the  cavity  boundary  and  satisfy  the  boundary  conditions  on  that 
approximate  surface.  However,  since  we  cannot  say  beforehand  that  our  approximate  location  of 
the  cavity  boundary  is  the  correct  one,  we  can  only  satisfy  the  dynamic  boundary  condition  (the 
condition  that  the  pressure  be  constant  on  the  boundary  and  equal  to  the  cavity  pressure)  in  the 
solution.  In  general,  the  kinematic  boundary  condition  (the  condition  that  the  normal  velocity 
vanish  on  the  cavity  surface)  will  not  be  satisfied.  However,  we  may  use  the  kinematic  boundary 
condition  to  update  the  location  of  the  cavity  surface. 

Since  the  first  guess  of  the  cavity  location  is  arbitrary,  it  is  numerically  convenient  to  first  solve 
the  boundary  conditions  on  the  wing  surface  and,  for  supercavities,  on  a zero-thickness  sheet 
downstream  of  the  wing.  To  find  the  nonlinear  solution,  one  simply  re-solves  the  problem  with  the 
new  body  shape  defined  from  the  first  iteration  (see  Figure  13).  The  fully  nonlinear  solution  is 
found  when  the  solution  converges  to  some  preset  tolerance.  In  previous  work  (Fine,  1992), 
extensive  studies  showed  that  the  accuracy  of  the  first  iteration  solution  relative  to  the  fully 
nonlinear  solution  was  very  good  over  a reasonably  large  range  of  operating  conditions.  As  a result 
(and  motivated  by  the  significant  computational  savings  of  solving  only  the  first  iteration)  LScav 
was  configured  to  compute  only  the  first  iteration  of  the  nonlinear  solution. 

Sample  LScav  results  are  shown  in  Figure  14,  which  indicates  the  convergence  of  the  cavity 
planform  for  a fin  with  a triangular  planform  and  a wedge-shaped  cross-section.  Table  1 shows  the 
convergence  of  the  corresponding  lift  and  drag  coefficients  and  the  cavity  volume.  Convergence  to 
three  significant  digits  occurs  with  about  1300  total  panels.  The  total  amount  of  memory  that  is 
required  for  this  size  problem  is  roughly  108  MB. 

LScav  has  been  validated  through  comparison  of  computed  lift  and  drag  coefficients  to 
measurements  made  in  a water  tunnel.  A sample  of  the  validation  is  shown  in  Figure  15. 

4 APPLICATION  OF  SLENDER-BODY-THEORY 

This  section  is  based  on  Varghese,  et  al,  (1997). 

In  slender-body  theory  (SBT),  the  slenderness  ratio,  e,  of  the  body  is  defined  as  the  ratio  of  the 
maximum  radius  to  the  total  length.  It  is  assumed  that  e « 1 and  that  the  cavitation  number  and 
the  inverse  square  of  the  Froude  number  are  of  order  e2ln  8 (Chou,  1974).  The  axial  gradient  of  the 
body-cavity  radius  should  be  of  the  same  order  as  the  slenderness  parameter,  restricting  the 
analysis  to  certain  classes  of  cavitators  using  this  approach. 

In  the  following  sections,  Chou’s  formulation  is  extended  to  subsonic  compressible  fluid  flows. 
Also  presented  are  selected  results  for  an  incompressible  fluid.  Chou’s  original  method 
incorporated  terms  accounting  for  an  axisymmetric  gravitational  field,  and  these  have  been  retained 
in  the  following  presentation.  Results  may  be  found  in  Varghese,  et  al,  (1997). 

4.1  Mathematical  Formulation 

The  present  SBT  analysis  assumes  an  axisymmetric  body  with  a conical  cavity  closure.  A salient 
cavity  detachment  point  is  specified  in  the  model  as  is  shown  in  Figure  16.  The  free-stream 
velocity,  is  at  zero  angle  of  attack.  The  body  length  is  the  cavity  length  is  £c,  and  the  total 
length  is  t.  In  the  following  formulation,  all  lengths  are  made  dimensionless  with  body  length, 
velocities  with  free  stream  velocity,  and  potential  with  body  length  and  free  stream  velocity. 

Incompressible  Flow 
The  total  potential  is 

^ = x + § , (4.i) 

where  <j)  is  the  perturbation  velocity  potential.  The  axial  and  radial  components  of  the  perturbation 
velocity  are 

u = — and  u = 

1 dx  r dr 


3(j) 


respectively. 


(4.2) 
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The  SBT  outer  expansion  results  in  Laplace’s  equation  for  the  perturbation  potential: 
<Kx,r)+—  + <!>„  = 0 . 

r (4.3) 

For  a continuous  distribution  of  sources  along  the  x-axis,  the  solution  for  equation  4.3  can  be 
written  as 


n i rc 

^=-i:l-f==7T 


4n  -J(x-tf+r2 


where  q(E,)  is  the  source  strength  defined  along  the  axis  at  location  E,.  From  the  SBT  inner 
expansion,  assuming  db/dx  = 0(e),  the  source  strength  is 

d(b2) 
q(x)  = 7l 

dx  (4.5) 

The  kinematic  boundary  condition, 

db 


is  formed  from  the  no-flux  condition  on  the  body  and  cavity  surfaces.  The  cavity  pressure  is 
assumed  to  be  constant.  The  dynamic  boundary  condition  is  derived  from  Bernoulli’s  equation: 

_ „ 2 2(*-l) 


cr  = -Cp  = 2 ux  + ur 


Here,  a is  the  cavitation  number,  Cp,  is  the  pressure  coefficient,  and 

F t— — 

b yfgt b 

is  the  Froude  number  based  on  body  length.  Here  g is  the  gravitational  acceleration. 

Let  £b  = b2h  and  C,c  = b2 . Substituting  expressions  for  the  radial  and  axial  velocity 

components  defined  at  the  body-cavity  surface,  in  the  dynamic  boundary  condition,  equation  4.7, 
results  in  a non-linear  integro-differential  equation: 


Wf  2(^1)  < 1<x<<, 
|(x-o!+cF  ldxl  F> 


where  bc  is  the  cavity  radius.  Integrating  by  parts  and  using  equation  4.5  for  the  source  strength 
yields  (Chou,  1974): 

d L df  —^-df  ^1 

1 f 1 r df  | dx1,  dx1,  , V 2(x  - I) 

2 \ J(x  - O'  -c  2 ,,  V(x  O1  1 C V(x  - ey  + c.  V*  + C ' dx ' K2 


db  \2  2(x  — 1) 


With  the  assumption  that  the  cavity  closes  as  a cone,  the  above  equation  can  be  written  for  a 
conical  body  as: 


- f d = 

2 , J(x  - O2  + C 


1 f dZ!  ^dxl  + 20f_l>l 

2 { J(  x - 0 ‘ + C < l\ 


(4.10) 
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In  this  equation,  the  right  hand  side  is  known  for  a given  cavity  shape.  In  the  current  approach,  the 
cavity  shape  is  approximated  and  updated  iteratively.  With  an  assumed  initial  cavity  shape, 
equation  4.10  is  solved  for  a new  cavity  shape.  For  that  procedure,  the  cavity  axial  length  is 
discretized  into  n panels.  In  each  panel,  a locally  quadratic  polynomial  is  used  to  approximate  C,c. 
Varghese,  et  al,  (1997)  presents  the  details  of  building  and  solving  a system  of  equeations  to  update 
the  discretized  cavity  geometry. 

Pressure  and  Drag  Coefficients  for  the  Case  of  an  Incompressible  Fluid 
Upon  convergence,  the  pressure  coefficient  along  the  body  is  computed  as  follows: 


Cp(x)  = - 


^(x-0  dH 

dC 

[(— tf+c„F 


2(x-l) 


0<x  <1. 


(4.11) 


The  pressure  coefficient  along  the  cavity  is  the  negative  of  the  cavitation  number.  The  base  drag 
coefficient,  based  on  the  maximum  body  cross-sectional  area,  is  found  by  integrating  the  pressure 
coefficients  over  the  body: 

■ [ C„  (. x ) — -dx+  c . 

l J p dr 

(4.12) 


Cdb  = 


In  Varghese,  et  al,  (1997),  viscous  drag  corrections  for  axisymmetric  bodies  were  developed 
based  on  three  different  methods:  (1)  Thwaites’  method  with  a linear  velocity  profile;  (2)  Thwaites’ 
method  with  a shear  function;  and,  (3)  the  method  of  Falkner  and  Skan.  Moderate  improvement  to 
the  comparison  of  the  predicted  total  cavitator  drag  to  experimental  data  was  obtained  for  certain 
cases. 

Subsonic  Compressible  Flow 

In  this  section  we  extend  the  previous  formulation  to  the  case  of  compressible  flows.  Results  may 
be  found  in  Varghese,  et  al,  (1997).  A method  is  proposed  in  Kirschner,  et  al,  (1995)  for  extending 
the  axisymmetric  boundary-element  method  presented  above  to  compressible  subsonic  flows. 

The  outer  expansion  for  the  perturbation  potential  for  slender  bodies  in  a compressible  fluid 
results  in  the  following  governing  equation  (Ashley  and  Landahl,  1965): 

(l-A4>ti+%  + 4>„.  =0  . 

r (4.13) 


where  the  hatted  variables  represent  the  physical  coordinate  for  the  subsonic  flows.  With  a 
transformation  to  x = x and  r = fir  where  (3  = -J  1 — , equation  4.13  reduces  to  Laplace’s 

equation.  Thus  the  solution  for  a continuous  source  distribution  is 


§(x,r) ■ 


r 


4n  °^(x-^f  +p2r2 


(4.14) 

Here  qix)  is  defined  by  equation  4.5.  The  radial  velocity  component  is  given  by  equation  4.6. 
However,  the  axial  velocity  component  depends  on  (3  as  follows: 


u 


X 


(*-^)2+/)2r2]2 


(4.15) 
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With  the  constant.cavity  pressure  assumption,  the  dynamic  boundary  condition  is  derived  from 
the  compressible  Bernoulli  equation: 


a = - 


yMi 


/-I 


(4.16) 


As  for  the  case  of  incompressible  flow,  substitution  of  equations  4.6  and  4.15  respectively  for  the 
radial  and  the  axial  velocity  components  of  equation  4.16  yields  an  integro-differential  equation: 


y-\ 


= i-(r-i)Mi 


1 f <t(Z)(x-Z) 


4 K I 


(x-|)2  + /32fc] 


if  dbc 
-rdg+ 

3 2 dx 


1<X<> 


(4.17) 


Integrating  by  parts  and  substituting  equation  4.5  for  the  source  strength,  the  governing  equation 
becomes 


d\ c 

c/(g2 


4 


^(x-zf+p2^  (r- l)M~ 

d2Zb 


Id 

r 

-l 

2 

" 

d%2 


dt;  d£c  | dQ b | 


dCb  | 


ofx 


c/.v 


J(x-t)2  + p2Zc  V-v2  + /*2fc 


(ff 


(4.18) 


For  a conical  body  with  a conical  cavity  closure,  this  equation  becomes 


d 2C 


^(x-d2+/32?c  (r-0^ 


7-1 

2 

7 

-1 

rf2?c 

A=2 


>J(x-Zt+P2(c  2?e 


(4.19) 


Pressure  and  Drag  Coefficients  for  the  Case  of  a Compressible  Fluid 
The  pressure  coefficient  along  the  body  is  given  by 


yMi 


cP  M= 


\-(y-\)Mi 


df b 
1 ,t  d% 
4Jo 


If  *b.T 

1 2{  dx  ) 


\x-^)2+p2Ch] 


7 

7-1 


-1,  0 <x  < 1. 


(4.20) 


Using  this  form  to  compute  the  pressure  coefficient,  the  base  drag  coefficient,  C^f, , for  the 
compressible  case  can  be  calculated  using  equation  4. 12. 


4.2  Selected  Results  - Slender  Body  Theory 

The  SBT  model  was  used  to  predict  supercavitating  flows  of  both  compressible  and  incompressible 
fluids  at  various  cavitation  numbers.  All  the  results  shown  here  are  for  a 10°-cone  angle.  The 
cavity  shapes  (scaled  to  show  cavities  enlarged  in  the  radial  direction)  for  cavitation  numbers 
ranging  from  0.005  to  0.05  are  shown  in  Figure  17  for  the  incompressible  case.  The  cavity  length 
increases  as  the  cavitation  number  decreases,  as  expected.  The  predicted  surface  pressure 
distributions  are  plotted  in  Figure  18  for  different  cavitation  numbers.  This  quantity  decreases  as 
the  cavitation  number  increases.  Results  have  been  compared  with  curve  fits  of  empirical  data 
(May,  1975).  In  Figure  19,  the  maximum  cavity  radius  has  been  plotted,  made  dimensionless  with 
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respect  to  the  radius  at  cavity  detachment.  Figure  19  illustrates  an  excellent  match  with 
experimental  values,  particularly  at  low  cavitation  numbers  where  SBT  is  more  valid. 

4.3  Summary  - SBT 

Formulations  have  been  presented  for  applying  Chou’s  SBT  model  to  both  incompressible  and 
compressible  flows.  Selected  results  have  been  presented  for  the  incompressible  case.  Some 
comparison  with  experimental  data  has  been  presented. 

In  Varghese,  et  al,  (1997),  additional  formulation,  results,  discussion,  and  comparison  are 
presented  for  the  effects  of  both  an  axisymmetric  gravity  field  and  viscous  drag. 

5 NUMERICAL  FLIGHT  AND  CONTROL  SYSTEM  SIMULATION 

In  this  section,  which  is  based  on  Kirschner,  et  al,  (2001),  we  describe  numerical  flight  simulation 
for  supercavitating  vehicles.  We  first  discuss  the  nonlinear  forces  acting  on  the  cavitator,  the  fins 
(if  present),  and  any  portions  of  the  hull  that  penetrate  the  cavity  boundary  during  excursions  from 
the  fully-enveloped  condition.  Next,  we  describe  simulation  of  vehicle  flight,  including  system 
stability  and  system  performance  during  execution  of  a banked  turn.  Without  control,  some  vehicle 
configurations  can  be  unstable,  whereas  a feedforward-feedback  strategy  can  control  some 
configurations  over  a range  of  turn  rates. 

5.1  Hydrodynamical  Model 

The  behavior  of  the  cavity  is  central  to  the  dynamics  of  a supercavitating  vehicle.  The  nominally 
steady  cavity  behavior  forms  the  basis  of  the  quasi-time-dependent  model  implemented  for  the 
current  investigation.  This  model  affects  not  only  the  forces  acting  at  the  nose  of  the  vehicle,  but 
also  the  immersion  of  the  fins  and  any  planing  forces  due  to  contact  between  the  body  and  the 
cavity  boundary.  It  is  the  cavity  that  makes  this  dynamical  system  not  only  nonlinear,  but  also 
dependent  on  strong  memory  effects  that  link  the  history  of  motion  to  the  instantaneous  cavity 
shape. 

The  three  parameters  that  characterize  time-averaged  or  quasi-steady  ventilated  flows  are  the 
cavitation  number,  the  cavity  Froude  number,  and  the  ventilation  coefficient.  The  cavitation 
number  is  defined  in  equation  2. 1 1 , and  has  been  discussed  extensively  in  the  formulations  above. 
The  cavity  Froude  number  (not  to  be  confused  with  the  related  Froude  number  based  on  body 
length  discussed  above)  and  the  ventilation  coefficient  are  defined  as 

A y 

(5.1) 

and 

c t Q 

CQ~ — T’ 

Vdl  (5.2) 

respectively.  Again,  g is  the  gravitational  acceleration.  The  cavitator  diameter  is  dc.  and  pc 
are  the  ambient  and  cavity  pressures,  respectively.  V is  the  vehicle  velocity  and  Q is  the 
volumetric  rate  at  which  ventilation  gas  is  supplied  to  the  cavity.  The  cavitation  number  expresses 
the  tendency  of  cavitation  to  occur  in  a flow,  and  is  thus  a principal  quantity  governing  cavity 
dimensions.  The  Froude  number  characterizes  the  importance  of  gravity  to  the  flow,  and  therefore 
governs  distortions  from  a nominally  axisymmetrical  shape.  The  ventilation  coefficient  governs 
the  time-dependent  behavior  of  the  cavity  as  ventilation  gas  is  entrained  by  the  flow. 

For  the  horizontal  trajectories  considered  in  these  lecture  notes,  the  cavitation  number  is  very 
low  (on  the  order  of  0.015  to  0.030)  and  the  cavity  Froude  number  is  typically  rather  high  (on  the 
order  of  50  to  100).  The  cavity  is  slender,  and  its  maximum  diameter  is  at  least  5 times  greater 
than  the  cavitator  diameter.  Cavity  tail-up  under  the  influence  of  gravity  is  small,  but  not 
insignificant  to  the  vehicle  design. 
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The  nomenclature  used  for  basic  vehicle  dynamics  during  the  remainder  of  these  lecture  notes 
is  presented  in  Figure  20.  The  x-axis  is  directed  forward  along  the  nominal  vehicle  long  axis,  y is 
to  starboard,  and  z completes  an  orthogonal  system  obeying  the  right-hand  rule.  The  origin  of  this 
coordinate  system  is  located  at  the  cavitator  pivot  point.  Following  Kiceniuk  (1954),  the  following 
relationships  have  been  employed  to  estimate  the  hydrodynamic  forces  and  moment  acting  on  the 
circular  disk  cavitator  selected  for  investigation: 

CD  (a,ac ) = = cD  o (1  + a)  cos2  ac 

1/2  pV2Ac 

CL  (cr , ac ) = 4C  — ac - = Cm (t  + a ) cos ac  sinac 
l/2  pV2Ac 


Cm  ) — 


Mc(o~,«c)  rQ 

l/2  pV2dcAc 


(5.3) 


where  Dc  and  Lc  are  the  drag  and  lift  forces  acting  on  the  cavitator,  Mc  is  the  hydrodynamic 
moment,  and  Ac  is  the  disk  area.  The  drag  coefficient  at  zero  angle  of  attack  and  cavitation 
number,  Coo , is  taken  as  approximately  0.805  based  on  empirical  data  (May,  1975). 

Various  analytical,  numerical,  and  semi-  and  fully-empirical  models  have  been  developed  that 
provide  estimates  of  the  maximum  cavity  diameter,  dm ax,  and  the  cavity  length,  tc . The 
analytical  formulae  of  Garabedian  (1956)  provide  useful  and  reasonably  accurate  approximations 
for  investigation  of  vehicle  dynamics: 

^max  _ \C D (<7,0) 

dc  V a (5.4) 


and, 

h.=  lcpMln  1 

dc  V cr  a 


(5.5) 


For  purposes  of  time-domain  simulation  of  vehicle  dynamics,  it  is  desirable  to  use  simple 
models.  In  a cylindrical -polar  system  with  coordinates  x,  r,  and  0 , the  nominal  axisymmetric 
shape  can  be  approximated  as: 

rc  _ f*  max 
d c 2 d c 


1- 


/ x/ dc  - i c I2dc 
y f c ^2 dc  j 


'2A 


where  rc{x,6 ) is  the  cavity  offset  (Munzer  and  Reichardt,  1950,  derived  from  a low-order  potential 
flow  model). 

The  effects  of  gravity  and  turning  on  this  generalized  ellipsoidal  cavity  shape  can  be 
approximated  as  a distortion  of  the  line  of  centers  (yc,zc)’  described  by 


yc(x)  = +- 


1 a 


( .. 


turn 


Kdc  j 


(5.7a) 


and 


:(*)=  — 


'jl'2 

Kdc  J 


(5.7b) 
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where  a mm  is  the  turn  acceleration  and  ag  is  the  tail-up  acceleration  of  the  cavity,  here 

approximated  using  the  semi-empirical  formula  recommended  by  Waid  (1957).  A binomial 
expansion  (retaining  the  second  term)  has  been  used  to  simplify  the  term  associated  with  turning, 
applying  the  principle  of  independent  development  of  cavity  sections  (Logvinovich,  1985)  and 
assuming  that  the  cavitator  traces  a circular  path  in  a horizontal  plane.  Although  the  distortion  of 
the  cavity  line  of  centers  due  to  turning  and  gravity  has  been  considered,  distortions  associated  with 
cavitator  lift  have  been  ignored.  Cavity  foreshortening  has  also  been  neglected,  and  distortions  of 
the  cavity  sections  from  their  nominal  circular  shapes  have  been  ignored  for  the  current 
investigation.  Under  these  assumptions,  the  cavity  locus  of  section  centers  is  coincident  with  the 
cavitator  trajectory. 

Equations  5.7  apply  to  a steady  horizontal  turn  with  gravity  acting  in  the  vertical  direction. 
Although  generalization  to  more  complicated  trajectories  would  not  be  difficult,  the  current 
formula  provides  insight  concerning  the  behavior  and  limitations  of  a supercavitating  vehicle  in  a 
turn.  The  consistent  formulation  shows  the  relative  influences  of  turning  and  gravity,  both  of 
which  become  increasingly  important  at  low  cavity  Froude  numbers.  To  the  order  of  terms 
retained  in  equations  5.7,  the  locus  of  cavity  section  centers  is  quadratic  in  the  axial  coordinate. 

The  importance  of  cavity  distortion  in  high  turn  rates  is  apparent  in  Figure  21,  which  presents 
results  for  an  extreme  turn  (in  this  case,  a 5-g  turn  that  is  probably  impractical,  but  is  illustrative  of 
the  cavity-body  interactions  important  to  the  problem).  The  nonlinear  behavior  is  clear  from  this 
figure.  Firstly,  as  the  cavity  distorts  from  a nominally  axisymmetric  configuration,  the  immersion 
of  the  fins  into  the  ambient  liquid  outside  the  cavity  becomes  asymmetric.  Therefore,  the  couple 
associated  with  symmetric  or  anti-symmetric  fin  trim  will  result  in  nonlinear  system  response. 
Secondly,  a system  designed  for  close  envelopment  of  the  vehicle  by  a nominally  axisymmetric 
cavity  (or  even  one  designed  for  cavity  tail-up  in  straight  and  level  flight)  will  be  subject  to 
nonlinear  forces  associated  with  afterbody  planing. 

Only  the  most  essential  aspects  of  the  unsteady  hydrodynamics  have  been  modeled.  Self- 
pulsation of  the  cavity  has  been  ignored,  but  changes  in  cavity  shape  due  to  cavitator  motion  have 
been  taken  into  account  via  a convolution  scheme.  The  fin  forces  are  evaluated  at  each  time  step 
via  a look-up  algorithm.  Thus,  the  three  basic  cavitating  fin  flow  regimes  (base-  and  partially 
cavitating  and  supercavitating)  contribute  to  the  unsteady  fin  forces,  but  various  higher-order 
hydrodynamic  effects  in  the  evolution  of  the  fin  cavity  are  not  modeled.  Most  importantly, 
hydrodynamic  added  mass  effects  have  been  ignored. 

Secondary  flows  - notably,  those  due  to  the  fins  and  to  afterbody  planing  - have  been  ignored 
in  the  current  effort,  although  the  theory  used  to  estimate  the  afterbody  planing  forces  accounts  for 
the  lowest-order  effect  of  the  spray  jet. 

In  the  nominal  condition  of  straight  and  level  flight,  the  supercavitating  vehicle  configurations 
of  interest  in  these  lecture  notes  are  supported  in  the  vertical  direction  by  a lift  force  acting  on  the 
cavitator  and  by  a combination  of  planing  and  lift  on  the  fins  near  the  after  end  of  the  vehicle.  The 
optimal  support  by  afterbody  planing  forces  has  not  yet  been  determined  with  certainty,  although 
the  results  presented  below  are  enlightening.  This  secondary  function  of  the  cavitator  as  a lifting 
surface  is  essential  to  operation  of  a supercavitating  vehicle.  Since  virtually  no  lift  is  provided  by 
the  gas-enveloped  hull,  the  vehicle  weight  must  be  supported  by  the  control  surfaces  both  forward 
and  aft  of  the  mass  center.  It  is  natural  to  provide  the  forward  lift  vector  by  orienting  the  cavitator 
at  an  appropriate  angle  of  attack.  Moreover,  accounting  for  changes  in  vehicle  mass  and  mass 
center  as  fuel  is  consumed,  and  maintaining  optimal  vehicle  orientation  with  respect  to  the  cavity 
during  a maneuver,  require  that  this  trim  angle  be  variable  and  controlled. 

The  forces  acting  on  the  fins  were  predicted  using  a fully  three-dimensional  boundary-element 
method  of  the  type  described  above,  supplemented  with  a viscous  drag  correction. 

Afterbody  forces  in  the  case  where  planing  on  the  cavity  surface  occurs  were  computed  using 
an  extension  of  Wagner  planing  theory  developed  by  Logvinovich  (for  example,  1980).  The  cavity 
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in  way  of  the  planing  region  is  approximated  as  a cylindrical  free  surface.  The  pressure  force 
normal  to  the  inclined  longitudinal  axis  of  the  cylindrical  hull  is  given  by 

f ' . <i\ 


Fb=n  p 


r2 


V2  sin 


ab  cos 


rb  + 2/?0 


1- 


h0  + A 


(5.8) 


where:  p,  is  the  hull  radius  (assumed  to  be  constant  over  the  planing  region);  rc , ab,  and  A are 
(respectively)  the  cavity  radius,  the  angle  of  attack  between  the  longitudinal  axes  of  the  body  and 
the  cavity,  and  the  difference  between  the  cavity  and  hull  radii  (all  averaged  along  the  planing 
region);  and  ho  is  the  immersion  depth  at  the  transom  measured  normal  to  the  cavity  centerline. 
Similarly,  the  moment  of  pressure  forces  about  the  transom  can  be  expressed  as 

M0=k  p 72  V 2 cos2  ab  — + — — — 

\b+2h0h0+A_  (5.9) 

A viscous  drag  correction  was  applied  for  the  current  investigation.  Forces  on  the  blast  tube  were 
neglected,  since  the  fluid  outside  the  time-averaged  cavity  boundary  in  that  region  is  characterized 
by  a significant  void  fraction,  so  that  the  overall  planing  forces  acting  on  the  reduced  sectional  area 
are  expected  to  be  rather  small. 


5.2  Candidate  Vehicle  Configuration 

Several  design  issues  must  be  addressed  to  completely  understand  the  dynamics  and  control  of 
supercavitating  vehicles.  One  concern  is  the  optimal  allocation  of  fin  lift  relative  to  planing  force 
to  support  the  vehicle  afterbody.  A related  question  is  the  need  for  fin  support:  if  stable  flight  can 
be  maintained  with  the  afterbody  completely  supported  by  the  planing  force,  the  elevators  and  their 
attendant  drag  and  complexity  become  unnecessary. 

Simulation  was  conducted  on  a vehicle  with  the  following  characteristics:  4.0  m in  length, 
0.2  m in  diameter,  and  with  a cavitator  diameter  of  0.07  m.  The  fins  were  located  3.5  m aft  of  the 
cavitator,  and  were  swept  back  1 5°.  Although  the  mass  properties  of  the  vehicle  will  change  as  the 
rocket  and  ventilation  fuels  are  consumed,  they  were  assumed  constant  for  purposes  of  the  current 
analysis.  A cruciform  fin  arrangement  has  been  assumed:  A pair  of  horizontal  elevators  provide 
some  component  of  steady  lift  to  support  the  afterbody  and  would  be  important  to  depth  changes, 
while  a pair  of  vertical  rudders  stabilize  the  vehicle  in  roll,  but  are  otherwise  deflected  only  during 
maneuvers. 

As  is  discussed  in  Kirschner  (2001),  under  certain  conditions,  a banked  maneuvering  strategy 
might  be  advantageous  for  a self-propelled  supercavitating  vehicle.  This  strategy  was  adopted  for 
the  investigation  described  herein. 


5.3  Simulation  of  Underwater  Flight 

The  current  investigation  focused  on  the  simplified  system  of  Figure  22,  for  which  high  frequency 
hydrodynamic  noise  excitation  and  other  noise  sources  have  been  ignored.  Moreover,  it  is  assumed 
that  knowledge  of  the  system  state  is  accurate,  eliminating  the  requirement  for  a state  estimator. 
Flowever,  the  nonlinearity  of  the  vehicle-cavity  system  has  been  retained.  Nonlinearity  enters  the 
system  via  the  vehicle,  fin,  and  cavity  behavior,  due  both  to  the  nonlinear  force  coefficients  and  to 
memory  effects  associated  with  cavity  evolution.  Control  has  been  implemented  via  a linear- 
quadratic  regulator  (LQR).  Cavity  memory  effects  and  the  discontinuities  in  the  force  coefficients 
require  incorporation  of  the  feedforward  model  (discussed  below)  to  maintain  acceptable  system 
performance. 

For  the  simplified  system,  the  input  to  the  regulator  is  the  error  vector,  x = - y , computed 

as  the  difference  between  the  commanded  state  vector,  jcd,  and  the  corresponding  actual  state 
vector,  y . Both  the  actual  and  commanded  state  vectors  have  the  form 


x y z (f>  0 iff 

surge  sway  heave  roll  pitch  yaw 


<—  displacements 


u v w p q r 
surge  sway  heave  roll  pitch  yaw 


rates 
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However,  several  of  these  states  might  not  be  measured  on  a vehicle.  Thus,  for  the  simulations 
discussed  below,  several  of  the  state  variables  were  eliminated  from  the  system.  Typically,  a 7- 
state  system  was  investigated,  assuming  that  both  rates  and  displacements  in  roll,  pitch,  and  yaw 
are  available  from  an  inertial  measurement  unit  and  that  some  measurement  of  vehicle  depth,  z , 
can  be  obtained. 

A classical  regulator  would  provide  a control  vector,  w , directly  to  the  hydrodynamic  control 
effectors  (in  this  case,  the  fins  and  the  cavitator): 

U\  U-2  W3  W4  W5 

starboard  lower  port  upper  cavitator 

However,  the  control  vector  of  interest  involves  both  the  feedback  output  from  the  LQR,  u b , and 
output  from  the  feedforward  model,  «f  : 

U = Hjj  + Uf 

(5.11) 


(5.10) 


Prior  to  flight  simulation  of  a supercavitating  body,  the  thrust  and  nominal  speed  are  specified. 
Values  of  the  state  and  control  vectors  are  then  determined  for  four  basic  flow  conditions: 

• Two  trim  conditions:  fin-supported  with  pitch  equal  to  zero,  and  planing-supported  with  fin 
deflection  (relative  to  the  vehicle  axis)  equal  to  zero;  and,  for  each  trim  condition, 

• Two  maneuvering  conditions:  straight  and  level  flight,  and  level  turning  at  a specified  constant 
turn  rate. 


Determination  of  these  equilibrium  conditions  involves  a nonlinear  root  search.  In  consideration  of 
the  discontinuity  that  can  occur  at  the  angle  where  planing  first  occurs,  a Nelder-Mead  scheme  has 
been  adopted. 

The  feedforward  model  investigated  in  this  effort  is  specialized  for  supercavitating  vehicle 
dynamics.  A scheme  has  been  selected  for  which  the  system  is  regulated  about  one  of  four  distinct 
fixed  points.  The  point  selected  depends  on  the  user-specified  operating  condition.  The  first  two 
fixed  points  are  defined  by  the  nominal  condition  of  straight-line  flight  for  either  fin  or  planing 
support  of  the  afterbody.  The  other  two  are  similar,  except  that  they  are  selected  as  equilibrium 
conditions  for  steady  turning  in  the  horizontal  plane.  Thus  the  output  from  the  feedforward  model  is 


U f G 


fin  supported,  zero  pitch,  straight  flight 
planing,  zero  fin  deflection,  straight  flight 
fin  supported,  banking  pitch,  turn  rate  c/tLim 
planing,  zero  fin  deflection,  turn  rate  uturn 


(5.12) 


and  feedback  control  is  optimized  to  regulate  motion  about  one  of  these  equilibrium  conditions  as 
the  vehicle  traverses  its  commanded  trajectory. 

The  nonlinear  equation  of  motion  involves  not  only  the  present  state,  but  also  the  history  of 
cavitator  motion,  c(t) : 

x = f(x,c(t-t\u) 

y = #(■*•) 

where  r represents  the  delay  between  motion  at  the  cavitator  and  action  along  the  afterbody.  For 
each  flow  condition,  s,  the  state-space  model  for  the  simplified  system  is  defined  in  the  usual 
manner  as: 


x = A?x+  B su 
y = C,x+  Dj« 


se  {l, — ,4} 


(5.13) 


The  Jacobian  matrices  are  computed  numerically  prior  to  flight  simulation  as  (for  example) 


Aij  = 


Ml 

dxj 


se  {h---, 4} 


(5.14) 


where  the  indices  i and  j vary  from  1 through  12  for  the  complete  state  vector  and  from  1 through 
5 for  the  control  vector  in  each  flow  condition.  For  the  simplified  system,  the  Jacobian  matrices  in 
the  second  of  equations  5.13  are  simply 

C,  = I .v  g {!,•••  ,4} 


(5.15) 
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(the  identity  matrix),  and 
D„  = 0 ss  {l, • • • ,4} 

(5.16) 

Regulation  for  the  general  LQR  is  measured  by  the  quadratic  performance  criterion 
j(u)=  j (jctQx  + 2xTNiq,  + u^Ru^Jdt 

t=o  , (5.17) 

where  the  user-specified  weighting  matrices,  Q , N , and  R , define  the  trade-off  between 
regulation  performance  and  control  effort.  For  the  current  investigation,  the  cross-weighting 
matrix  N was  taken  as  zero,  as  were  all  off-diagonal  elements  of  the  state  and  control  weighting 
matrices,  Q and  R.  A MATLABlil  Control  System  Toolbox  routine  is  used  to  determine  the  gain 
matrices,  Ks,  that  minimize  the  cost  function,  J[u ) , for  use  in  the  state  feedforward- feedback 
control  law  (MATH  WORKS,  1998)  for  each  of  the  4 fixed  points  considered.  The  control  law  is 
then  expressed  as: 

U = Mfs  -K„x  S E -fl,  — ,4} 

1 . (5.18) 

A substantial  effort  was  devoted  to  testing  various  combinations  of  the  weighting  matrices  Q and 
R to  achieve  acceptable  performance,  although  a more  systematic  and  rigorous  effort  is  warranted. 
The  weighting  and  gain  matrices  are  all  computed  prior  to  flight  simulation. 

A 4th-order  Runge-Kutta  scheme  was  used  to  integrate  the  equations  of  motion  in  time.  At 
each  time  step  the  cavity  geometry  is  updated  to  account  for  the  motion  of  the  cavitator.  The  forces 
acting  on  the  cavitator,  fins,  and  afterbody  are  then  updated,  and  vehicle  accelerations  in  each  of 
the  6 degrees  of  freedom  are  computed.  The  state  and  error  vectors  are  updated.  The  control 
vector  is  then  updated  using  the  gain  matrix  appropriate  to  the  specified  desired  operational 
condition  and  the  control  equation  is  applied. 

5.4  Selected  Results  - Simulations 

For  all  cases  discussed  below,  the  operating  speed  and  cavitation  number  were  75  m/s  and  0.03, 
respectively. 

Figure  23  shows  the  total  lift  margin  (the  difference  between  the  total  lift  and  the  vehicle 
weight)  as  a function  of  vehicle  pitch  and  elevator  deflection  angle.  The  cavitator  angle  for  each 
curve  was  selected  to  provide  equilibrium  for  a pitch  angle  of  zero.  Equilibrium  is  achieved  in  the 
fin-supported  case  at  elevator  deflection  angles  between  approximately  0.30°  and  0.55°;  below 
elevator  deflection  angles  or  approximately  0.30°  the  afterbody  is  primarily  supported  by  planing 
forces.  Note  the  discontinuity  in  lift  curve  slope  when  afterbody  planing  first  occurs  (at  a pitch 
angle  of  approximately  0.23°);  this  represents  a tremendous  increase  in  the  effective  system  stiffness 
and  the  characteristic  response  frequencies.  The  small  fin  deflections  required  to  produce  a 
substantial  change  in  the  forces  acting  on  the  body  might  present  an  additional  challenge  for  control. 

The  open-  and  closed-loop  eigenvalues  of  the  system,  linearized  about  the  fixed  points 
associated  with  the  fin  and  planing  supported  cases  for  straight-line  flight,  are  shown  in  Figure  24. 
For  both  types  of  afterbody  support,  the  uncontrolled  system  has  one  unstable  mode.  Preliminary 
analysis  suggests  that  this  is  a corkscrew  motion  associated  with  cavity  tail-up  under  the  influence 
of  gravity:  roll  and  sway  are  unstably  coupled  due  to  asymmetry  in  the  fin  immersion  as  the  body 
rolls  about  its  long  axis,  which  is  eccentrically  located  with  respect  to  the  line  of  cavity  centers  in 
way  of  the  fins.  This  effect  can  be  eliminated  in  the  controlled  system,  provided  some 
measurement  of  depth  is  available;  if  not,  the  unstable  mode  is  present  even  for  the  controlled 
system.  As  expected,  the  most  important  remaining  difference  between  the  controlled  and 
uncontrolled  systems  is  a significant  increase  in  the  damping  of  each  mode  in  the  controlled  case. 
Similarly,  the  primary  difference  between  the  fin  and  planing  supported  cases  is  a significant 
increase  in  the  frequency  of  the  oscillatory  modes  in  the  planing  supported  case.  In  all  cases, 
several  lightly  damped  modes  are  present. 

The  vehicle  motion  predicted  by  time-domain  simulation  of  free  decay  in  pitch  for  the  fin- 
supported  case  in  straight  and  level  flight  is  presented  in  Figure  25.  Two  types  of  behavior  are 
presented:  essentially  linear  behavior  associated  with  a small  initial  pitch  rate  perturbation  of  l°/s, 
and  the  clearly  nonlinear  behavior  due  to  a large  perturbation  of  20°/s.  In  both  cases,  control  is 
required  to  recover  from  the  perturbation  and  return  to  the  original  depth.  More  interestingly,  the 
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linear  behavior  is  non-oscillatory  for  both  the  controlled  and  uncontrolled  cases,  whereas,  in  the 
nonlinear  case,  the  perturbation  initiates  a rather  violent  pitch  oscillation  at  a frequency  of 
approximately  9.5  Hz  that  is  modulated  by  a second,  gentler,  higher-frequency  oscillation.  The 
shapes  of  the  peaks  and  troughs  indicate  that  the  vehicle  is  skipping  off  the  bottom  of  the  cavity  in 
the  nonlinear  case. 

Simulation  of  vehicle  motion  in  steady,  planing-supported,  straight  and  level  flight  predicts 
behavior  such  as  that  shown  in  Figure  26.  It  can  be  seen  that  pitch  increases  with  time,  but  that  a 7- 
state  controller  both  eliminates  the  tendency  of  the  vehicle  to  dive  and  reduces  the  violent  tail-slap 
behavior.  The  frequencies  characterizing  the  two  types  of  response  are  markedly  different:  the 
uncontrolled  case  is  dominated  by  a relatively  long-period  skipping  mode,  whereas  control 
eliminates  this  mode,  leaving  the  higher  modulating  frequency.  The  planing  force  history  indicates 
that  skipping  occurs  in  both  cases  - apparently  the  controlled  response  during  the  planing  portion 
of  the  cycle  is  sufficient  to  mitigate  the  unsteady  cavity  behavior  in  such  a manner  that  the 
afterbody  excursions  are  not  so  extreme. 

5.5  Summary  - Simulation 

This  section  presented  strategies  for  the  control  of  the  highly  coupled  nonlinear  system 
comprising  a supercavitating  vehicle.  A simple  hydrodynamic  model  was  implemented  for 
simulating  the  behavior  of  such  a system.  In  consideration  of  the  nature  of  the  forces  acting  on  a 
gas-enveloped  body,  and  in  order  to  maintain  mechanical  simplicity  of  the  cavitator,  a banked 
maneuvering  strategy  was  investigated.  Results  of  dynamical  simulation  for  a specific  vehicle 
were  presented  both  for  uncontrolled  flight  and  with  LQR-based  feedforward-feedback  control. 

The  system  behavior  is  dominated  by  the  distinct  change  in  the  nature  of  the  forces  as  the 
afterbody  moves  between  a planing  and  a non-planing  condition.  The  system  eigenvalues  strongly 
depend  on  the  type  of  afterbody  support  specified  for  vehicle  operation,  but  that  (under  the 
assumptions  made  for  the  simplified  hydrodynamics  model)  a 7-state  controller  can  eliminate  the 
most  undesirable  behavior  in  either  case.  Additional  discussion,  including  description  of  the 
simulation  of  banked  maneuvers  in  steady  turns,  may  be  found  in  Kirschner,  et  al,  (2001). 
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Figure  1.  Two  classes  of  supercavitating  high-speed  bodies:  (a)  Free-flying  projectiles;  (b)  Self-propelled 

vehicles 
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Figure  6.  Body-cavity  profiles  for  a 15°  cone  for  dimensionless  cavity  lengths  shown  (scaled  for  viewing). 
The  corresponding  predicted  cavitation  numbers  (decreasing  with  length)  are  0.501,  0.070,  0.038,  and 

0.017 
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Figure  7.  Predicted  drag  coefficients  versus  cavitation  number  for  a 15°  cone  at  dimensionless  cavity 
lengths  shown.  The  empirical  curve  fit  of  May  (1975)  is  based  on  results  of  experiments  by  various 

researchers. 
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Figure  8.  Predicted  cavity  diameter  versus  cavitation  number  for  a 15°  cone  at  dimensionless  cavity 
lengths  shown.  The  empirical  fit  of  May  (1975)  is  based  on  results  of  experiments  by  various  researchers. 
Solid  markers  are  based  on  the  computed  drag  coefficient;  open  markers  are  based  on  the  drag  coefficient 
given  by  the  empirical  curve  fit  of  May  (1975).  The  function  of  diameter  is  given  in  the  text  at 

equation  3.18. 
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Figure  10.  The  initial  boundary-value  problem  for  unsteady  boundary  element  simulations. 
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Figure  11.  Normal  velocity  at  start-up,  with  a homogenous  initial  condition  for  the  cavity  potential  and 

an  assumed  initial  geometry. 


Figure  12.  The  three  regimes  of  cavitating  wedges,  fully-,  partially-  and  base  cavitating. 
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Figure  13.  Sketch  demonstrating  the  difference  between  the  first  iteration  and  nonlinear  solutions. 

LScav  solves  for  the  first  iteration  solution. 
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Figure  14.  Convergence  of  cavity  planform  a fin  with  45°  sweephack  and  7.5°  wedge  angle  (note:  the  plot 

is  not  drawn  to  scale). 


Table  1.  Sample  LScav  results  showing  convergence  of  cavity  volume,  lift,  and  drag  coefficients  with 

number  of  panels 
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Figure  15.  Comparison  of  measured  and  computed  lift  and  drag  coefficients  (measured  data  provided  by 

D.  Stinebring,  ARL/PSU) 
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Figure  16.  Physical  problem  addressed  by  slender  body  theory. 
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Figure  1 7.  Cavity  shapes  for  different  cavitation  numbers  for  axisymmetric  supercavitating  flow  past  a 

10°  cone  in  an  incompressible  fluid. 
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Figure  18.  Predicted  Surface  Pressure  Distributions  for  Various  Cavitation  Numbers  for  Axisym metric 
Supercavitating  Flow  Past  a 10°  Cone  in  an  Incompressible  Fluid 


Figure  19.  Cavity  Maximum  Radius  Versus  Cavitation  Number  for  Axisymmetric  Supercavitating  Flow 

Past  a 10°  Cone  in  an  Incompressible  Fluid 
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Figure  21.  Cavity  behavior  in  an  extreme  turn 
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Figure  22.  Simplified  dynamical  systems  model  of  a supercavitating  vehicle. 
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PITCH  (°) 

Figure  23.  Lift  margin  as  a function  of  pitch  and  elevator  angles. 


Figure  24.  Open-  and  closed-loop  eigenvalues  for  two  modes  of  afterbody  support. 
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Initial  Perturbation  in  Pitch  Rate  = 1 °/sec 
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Figure  25.  Fin-supported  straight  and  level  flight  with  small  and  targe  perturbations. 


